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Abstract Glioblastoma is a rapidly evolving high-grade astrocytoma that is 
distinguished pathologically from lower grade gliomas by the presence of necro- 
sis and microvascular hipcrplasia. Necrotic areas arc typically surrounded by 
hypercellular regions known as pseudopalisades originated by local tumor ves- 
sel occlusions that induce collective cellular migration events. This leads to the 
formation of waves of tumor cells actively migrating away from central hypoxia. 
We present a mathematical model that incorporates the interplay among two 
tumor cell phenotypes, a necrotic core and the oxygen distribution. Our simu- 
lations reveal the formation of a traveling wave of tiimor cells that reprodiiccs 
the observed histologic patterns of pseudopalisades. Additional simulations of 
the model equations show that preventing the collapse of tumor microvessels 
leads to slower glioma invasion, a fact that might be exploited for therapeutic 
purposes. 
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1 Introduction 



Malignant gliomas are the most common and lethal type of primary brain 
tumor. Survival for patients with glioblastoma multiforme (GBM), the most 
aggressive and prevalent WHO grade IV astrocytic glioma ( [Louis et al.||2007 



Wen and Kesari 2008 1, although individually variable, is about 12 months after 



diagnosis, using the standard of care which includes surgery to resect as much 
tumoral tissue as possible, radiotherapy and chemotherapy (temozolamide) 



(Huse and Holland 2010 Van Meir et al. 2010 Clarke et al. 2010 Reardon et 



al. 2011 Pruitt et al. 2011 1. GBM is a rapidly evolving astrocytoma that is dis- 



tinguished pathologically from lower grade gliomas by the presence of necrosis 



in the central regions of the tumor and microvascular hyperplasia ( Louis et al. 



2007). Another distinctive signature of GBM is the fact that necrotic foci are 



frequently surrounded by hypercellular regions known as pseudopalisades. An 
example of these hypercellular perinecrotic structures is shown in Figure [ij 

The fact that GBM is one of the most highly vascularized human tumors 
seems to contradict that its microcirculation is functionally very inefficient 
( Jensen||2009 ) . However, it is well known that the tumor microenvironment is 
remarkably different from that of normal tissue; its chaotic vasculature leads 
both to a decrease in the supply of oxygen and essential nutrients as well as 
in the removal of waste products. 

Tumor hypoxia (deficit in oxygenation) is a feature encountered in most 
solid tumors ^Semenza et~al]|2003[ |Vaupel|[2008l [Bristow and Hill||2008[ [Pe^ 
2008), albeit with an incidence and severity that varies among 



whirst et al. 



patients. It is generally recognized as a negative clinical prognostic and predic- 
tive factor owing to its involvement in various cancer hallmarks such as resis- 
tance to cell death, angiogenesis, invasiveness, metastasis, altered metabolism 
and genomic instability ( Hanahan and Weinberg|2011 Wilson and Hay|2011 1. 
Hypoxia displays a central role in tumor progression and resistance to ther- 



apy (chemo- and radioresistance), specially in GBM ( Jensen[2009 Pope et al. 



2011). Low oxygen levels are observed both in the near vicinity of the tumor 



vasculature, as acute or cycling hypoxia (Bristow and Hill 2008 Dewhirst et 



al. 2008 1, or at longer distances of about 150 fim from the feeding blood vessels 



giving rise to areas of chronic or diffusion-limited hypoxia (Bristow and Hill 



2008 ) . Angiogenesis emerges then in response to the imbalance of proangio- 



genic growth factors that are released by hypoxic cells in the tumor such as 
vascular endothelial growth factor (VEGF) relative to anti-angiogenic growth 



factors (e.g., angiostatin) present in the tumor microenvironment (Carmeliet 
and Jain|20TI||Ebos and Kebel|20l"T| ). The end result of VEGF signaling in tu 



mors is the production of immature, highly permeable blood vessels with sub- 
sequent poor maintenance of the blood brain barrier and parenchymal edema 
( Jain et al.|2007 ). Thus, microvascular hyperplasia may coexist with moderate 
or high levels of hypoxia. Actually, for GBM it has been shown that microvas- 
cular hyperplasia is spatially and temporally associated with pseudopalisading 
necrosis and is believed to be driven by hypoxia-induced expression of proan- 
giogenic cytokines such as VEGF ( Rong et al.||2006" ). 
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Fig. 1 Pseudopalisading structures in Glioblastoma Multiforme. This liistopatlio- 
logic sample shows the typical features of pscudopalisades in GBM. Hypercellular regions 
appear around centrally degenerating vessels (star). The standard interpretation ( |Brat et] 
|al.|2004"] | is that cells migrate away from the degenerated vessels (hypoxic) areas and move 
towards functional vascular supplies such as the one marked with an arrow head in the 
photograph. 



In general, GBMs are pathologically heterogeneous with central regions of 
necrosis surrounded by pseudopalisading cells and hypervascularized regions 
under moderate levels of hypoxic stress ( Gorin et al. |2004, ) . In addition, the 
periphery of the GBM typically shows tumor cells infiltrating into the nor- 
mal brain. This invasion often extends to the adjacent cortex and through 
the corpus callosum into the contralateral hemisphere, being a major reason 
for treatment failure: migrating cells are not eliminated by surgical resection 
and cause tumor recurrence ( Onishi et al.pOlT ). Thus, high grade astrocytic 
tumors typically recur in less than six months after surgery ( |Giese et al.[2003| . 

Emerging evidence from clinical observations suggests that the origin of 
pseudopalisade formation is a vaso-occlusive event where the local blood ves- 
sels no longer provide the necessary oxygen supply to the tumor microenviron- 
ment, resulting in a population of migrating hypoxic cells responsible for the 
formation of the pscudopalisades (see Figure [T]) . Establishing a clear mecha- 
nism underlying pseudopalisade formation and its role in necrosis is critical 
since their formation could potentially precede all forms of coagulative necrosis 
in GBM (? Brat and Van Meir|20"04{|Brat et al.|2004[ ), responsible for the broad 
necrotic areas detected, for example, by standard magnetic resonance imaging 
of brain tumors. In fact, pathologic examinations indicate that hypoxia and 
necrosis development within high grade astrocytomas could arise secondary 
to vaso-occlusion. Pscudopalisades would then represent a wave of tumor cells 
migrating away from central hypoxia ( Rong et al.||2006 |. 

Experimental studies have hypothesized that pseudopalisade formation 
could provide a mechanism accounting for the rapid clinical progression of 
GBM, acting as a link among the underlying vascular damage, the develop- 
ment of hypoxia and hypoxia-induced angiogenesis, leading to necrosis and 
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Fig. 2 Glioblastoma microenvironment modeled in this work. The system modeled 
consists of two blood vessels and an evolving embedded population of tumor cells wliose phe- 
notype changes according to the levels of oxygenation. Upon thrombosis of one of the blood 
vessels (the right one), due to excessive tumor cell growth and the secretion of thrombotic 
factors, decreasing levels of local oxygenation arise that promote a massive migration of 
hypoxic glioma cells (labeled as II) towards better oxygenated areas. In the process, the hy- 
poxic glioma cells infiltrate among the normoxic glioma cells (labeled as I) , which are closer 
to a functional vessel (the left one), creating a transient hypercellularity region (a pseu- 
dopalisade). As the pseudopalisading front of glioma cells enlarges around the thrombosed 
vessel, perivascular necrosis becomes more prominent. 



accelerated outward tumor expansion in the form of a migrating wave of cells. 
Previous works have proposed that pseudopalisade formation would be the 
result of a multistep process ( Brat et al. 2004[ |Rong et"al 



2006 Brat and 



Van Meir||2004 ) : First, in anaplasic astrocytoma (WHO grade III) or de 



GBM in its initial stages, tumor cells proliferate and infiltrate through the 
parenchyma and receive oxygen and nutrients via the intact native blood ves- 
sels. Secondly, vascular insult occurs as a result of uncontrolled tumor growth 
causing endothelial injury and vascular leakiness. Both endothelial injury and 
the expression of procoagulant factors by the neoplasm result in intravascu- 
lar thrombosis and increasing hypoxia in the regions surrounding the vessel 
(Rahman et al. 2010). Subsequently, tumor cells begin to migrate away from 
hypoxia, creating a peripherally moving wave that is seen microscopically as 
pseudopalisading cells. This leads to an expansion of the zone of hypoxia and 
central necrosis, whereas the hypoxic tumor cells in pseudopalisades secrete 
proangiogenic factors (e.g. cytokines VEGF and IL-8). 

In this work we study the formation of hyper cellular regions in perinecrotic 
areas in high grade gliomas. We put forward a mathematical model that incor- 
porates the spatio-temporal interplay among two tumor cell phenotypcs corre- 
sponding to well oxygenated (normoxic) and hypoxic cells, a necrotic core and 
the oxygen distribution. Our physio-pathologic scenario is depicted in Figure 
[2] and consists of a tumor cell population embedded within two blood vessels 
arranged in a flat domain. Our numerical simulations reveal the formation of 
a superimposed traveling wave of hypoxic cells that qualitatively reproduces 
the experimentally observed patterns and may be used to obtain some hints 
for the timescales of palisade formation, lifetime and persistence among other 
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prognostic metrics. We also explore the dynamics of tumor spreading under 
delayed vascular injury and show that somehow, contrary to naive intuition, in 
the framework of our model preventing vessels from breaking leads to slower 
tumor invasion speeds that might be explored for therapeutic purposes. 



2 The Model 



2.1 Glioblastoma compartment dynamics 



Our modeled system comprises three different compartments; two different 
coupled tumor cell subpopulations corresponding to the two dominant phe- 



notypes well described in GBM ( DeBerardinis et al. 2007 Giese et al. 2003 
Keunen et al.|20lT Onishi et al. 2011 1 and a necrotic compartment. Of the two 



populations, the first one, C„, consists of normoxic, proliferative, less mobile 
cells (typically located close to functional blood vessels), whereas the second 
one, C/j, is composed of hypoxic, less proliferative and more mobile cells. Both 
phenotypes compete for space and resources (in our case oxygen). The cell 
loss is included in a third compartment, Cd, of necrotic tissue. The equations 
governing the interplay among these three densities are as follows 



dt 



dt 



DnV^Cn + — f 1 (jiM) I ^ ShnCh 



1 

Tnh 



1 

Th 



1 - 



Cn + Ch + Cd 



(J(M) 



Thn 



Thn 



H SnhC'n ShdCh-, 

T nh T hd 

TTT = ShdCh- 

at Thd 



dC, 



(la) 

(lb) 
(Ic) 



The first terms in Eqs. (la) and (lb) account for the cellular motility. Migra- 



tion in gliomas is not simple and in fact many have proposed that the highly 
infiltrative nature of human gliomas recapitulates the migratory behavior of 
glial progenitors ( Suzuki et al.||2002 Dirks||2001 1. Here we assume, as in most 



models ([lYieboes et al.||2007[ [Swanson et ai.|2008[ [Bondiau et al.|2008[|Eiken 



berry et al.||2009| [Konukoglu et al.||2010[ |Rockne et al.||2010| |Perez-Garcia et 
al.||2011 ) , that glioma cell invasion throughout the brain is basically governed 
by a standard Fickian diffusion process. 

It is believed that glioma cells -as many other types of cells- do not move 
and proliferate simultaneously ( Giese et al.|2003 ). This is referred to as the mi- 
gration/proliferation dichotomy and leads to the fact that highly motile cells 
should exhibit low proliferation rates. The switch between proliferative and 



invasive phenotype cannot be only mutation driven (Hatzikorou et al. 2010 



Onishi et al. 2011 ) and it has been suggested that invasive glioma cells are 
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able to revert to a proliferative cellular program and vice versa, depending on 
the environmental stimuli (Giese et al. 2003 Keunen et al. 2011). Recently, 



it has been proposed that oxygen concentration may be one of those stim- 
uli able to drive this transformation and thus for each oxygen level, there 
exists a dominant (fittest) tumor cell phenotype that corresponds to certain 
ratio proliferation/migration rates (Giese et al. 2003). Several studies have 
linked hypoxia to the invasive behavior of different types of tumors and their 



relationship with metastasis and negative prognosis (Bristow and Hill 2008 



Kalliomaki et al. 2009 Elstner et al. 2007). Mathematically, the observation 



that the hypoxic 


phenotype is more migratory than the normoxic one ( . 


Berens 


and Giese 


1999 


Giese et al. 


2003 


Bristow and Hill 


2008 


Gorin et al. 


2004) 



translates to choosing the hypoxic cell diffusion coefficient Dh to be larger 
than the normoxic one I?„. 

The second terms in Eqs. (la) and (lb) employ a classical logistic growth 
for the tumor cell populations with proliferation times t„ and t/j, respectively. 
Thus, the net proliferation is lower in regions of high cell density where C„ -f 
Ch + Cd is close to the maximum density capacity C*-*"^-*, than in regions of 
low cell density where Cn -\- Ch + Cd is much smaller than C'*^-*. As to the 
parameters t„ and t/,,, the migration-proliferation dichotomy suggests that 
Tn < Th- Since growth is assumed to be space- limited we incorporate also the 
necrotic tissue density into the saturation terms. 

The third terms in Eqs. ( la ) and ( lb ) are essential ingredients of our model 
since they represent the phenotypic switch under the effect of changing mi- 
croenvironmental conditions. The switch functions Snh, Shn depend on the 
oxygen concentration (to be specified below) and account for the fact that 
under low oxygen conditions normoxic cells change their phenotype to the hy- 
poxic one with a rate 1 / t„^ , whereas when exposed to sufficiently high oxygen 
concentrations hypoxic cells recover their oxic phenotype with a rate 1 / t/i„ . 

Equations ([T]) also include hypoxic cell loss by death due to anoxia [cf. 
last term in Eq. ( |lb[ )] and how these lost cells feed the necrotic tissue [cf. Eq. 
(Icl]. The last term in Eq. (lb) represents (anoxic) cells dying when the oxy- 
genation conditions are adverse enough as measured by the switch function 
Shdi with a rate l/rhd- It should be emphasized that cell loss in glioma occurs 
both by apoptosis (representing death by controlled suicide) and by necro- 
sis (death by uncontrolled damage). Necrosis has historically been considered 
to be an inappropriate or accidental death that arises under conditions that 
are extremely unfavorable such as those incompatible with a critical normal 
physiological process. More recently, a number of studies have suggested that 
necrosis is also a regulated process that can be modulated. For example, in- 
duction of necrosis seems to be dependent on cellular energy stores, such as 
NAD and ATP. Furthermore, cell stress and cell signaling including oxidative 
stress, calcium levels and p53 activation have been shown to influence lyso- 
somal membrane permeability and survival ( Boya and Kroemer||2008 ) . When 
apoptosis sets in, the cell membrane collapses and the cell shrinks, while when 
necrosis occurs the cell keeps its shape and thus still occupies a physical space 
( Hotchkiss et al.||2009 ). This leads to the fact that cell loss due to apoptosis. 
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not leading to the fill up of space, can be simply incorporated into effective 
proliferation cell rates in Eqs. (la) and (lb). However, cell death by necrosis 



happens at a different rate and, in addition, results in occupied space; this 
is the reason for incorporating the population Cd into the proliferation lim- 
iting terms in Eqs. (la) and (lb). Overall, levels of apoptosis are generally 



low in malignant gliomas; apoptosis accounting for only a slim minority of cell 



death. Furthermore, apoptotic rates do not correlate with prognosis (Migheli 
et al.|1994 SchifFer et al.|1995 ). In contrast, coagulative necrosis represents the 



majority of cell death in GBMs and its degree is inversely related to patient 



survival, as several shown by Nelson et al. (1983) and Lacroix et al. (2001). 



2.2 Microenvironment oxygenation 



Though the GBM microenvironment is highly heterogeneous at the cellular 



and molecular levels (Bonavia et al. 2011), two of the main chemical agents 
implicated in its growth and metabolism are oxygen and nutrients, mainly 



glucose and lactate (Mendoza et al. 2011 Griguer et al.||2005 2008 Seyfried 
et al.|20lT ) . A more comprehensive description of the GBM microenvironment 
should, a priori, take into account the interplay not only of oxygen but of nu- 
trients as well, at the expense of increasing the complexity of the mathematical 
model. Here we show that it suffices to take oxygen as the key chemical agent 
driving the collective cell migration dynamics to understand pseudopalisade 
formation in GBM. Moreover, glucose is typically less scarce than oxygen even 
at long distances from the blood vessels and can be replaced by other fuels 



such as lactate or, in a smaller fraction, by pyruvate, malate or glutamine (De- 



Berardinis et al.|2007[ [Mendoza et~aL][20TT] [Beckner et "ar][2005l [Grillon et a~ 



2011 ). For instance, sequestration of lactate and pyruvate in the extracellular 



matrix has been observed in C6 gliomas in a rat brain by Grillon et al. (2011 ) 



and correlates with the persistent acidity of the GBM microenvironment. 

Oxygen tissue heterogeneities are recognized to be a very relevant factor in 
gliomas and many other tumor types ( Evans et al.|2004 ). The spatio-temporal 
variation of oxygen distribution, driving the presence of various phenotypic 
tumor subpopulations, would help to explain the diversity of responses ob- 
tained from the same treatment employed in patients with apparently similar 
tumors but irrigated in a different way. Here, a Michaelis-Menten type kinetics 
is used for the uptake of oxygen (Patel et al. 2001 Ferreira et al. 2002), but 



now encompassing the feedback by the normoxic and hypoxic cells 



do 



dt 



O2- 



OLnCn + OlhCh 



o 



(T) 



0-2 



(2) 



The first term in Eq. ^ accounts for the oxygen diffusion in the brain tissue 
assuming a homogenous and isotropic diffusion coefficient D02 for simplic- 
ity. Although oxygen passes successively through the intracellular fluid, cell 
membranes and cytoplasm, all having abrupt spatial variations, employing an 
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Oxygen 




+ proliferation - proliferation 

- diffusion + diffusion 



Fig. 3 Oxygen concentration influences the phenotype of a tumor ceU. Depending 
on the oxygen conecntration various switching ineelianisms arise coupling the populations: 
Normoxic to hypoxic, S„fi> hypoxic to normoxic S^n a-nd hypoxic to necrotic Sfni- High 
oxygen levels favor the existence of more proliferative phenotypes which are less mobile. 
On the contrary, cells respond to low oxygen concentrations by expressing less proliferative 
phenotypes which are more motile, and finally, hypoxic cells experiencing persistent anoxia 
will eventually die ( |Brat et al.|2004^ . 



average diffusion coefficient for oxygen lias been proven to be a good approxi- 
mation in previous works ( Tannock|1972 Pogue et al.|2001 ). Tfie second term 
models the oxygen consumption by normoxic and hypoxic cells at rates a„ and 

(T) 

ah, respectively. The saturation Michaelis-Menten constant corresponds 
to the oxygen concentration level at which the reaction rate is halved. 

When modeling the oxygenation dynamics, we will assume that normoxic 

cells switch to the hypoxic phenotype with a characteristic time t„^ when the 

(s) 

local oxygen concentration is below a critical value O2 , such that hypoxia- 
inducible factor la (HIF-la) is not degraded and the glycolytic and angiogen- 



esis mechanisms initiate (Jewell and Gassmann 20011. The switching of cells 



from the hypoxic to the normoxic state Shm as described by Eqs. (la) and 



(lb), occurs with a characteristic time r^„ and is a reversible process Snh] if 



the local oxygen sensed by a hypoxic cell is above 02^'' , the cell will undergo a 
phenotypic switch to the normoxic state (as it satisfies O2 > O2 ) and vicev- 
ersa. Finally, hypoxic cells experiencing persistent anoxia will eventually die 
when the oxygen concentration is below a threshold value 02^^ with a charac- 
teristic time Thd- The specific form of the three step-like switch functions are 



Shni02) = 
Snh{02) = 
SUO2) - 



1 + tanh 



1 — tanh 



O2-O2 
AO2 



{sy 



l-tanh ^ 
I AO2 



{sy 



AO2 



(D)' 



(3a) 
(3b) 
(3c) 
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The extra parameter AO2 accounts for the fact that transitions occur around 
the critical values Oj^' and 02^"* with a characteristic range (or width). In our 
calculations we have taken the same small width AO2 for all of the transitions. 
The change from the hypoxic state to necrosis is irreversible and takes place 
in an interval centered around the size of this transition domain being 

proportional to A02- 

Although it is possible to incorporate the oxygen supply by resorting to 
extra source functions, here we will use a more direct approach and oxygen 
will be provided initially through the two capillaries located at the boundaries 
of the tumor domain (see Figure [2]) . As soon as one of the two blood vessels 
becomes impaired by a thrombus, the remaining functional one will be the 
sole source of oxygen. 



2.3 Parameter estimation 



We resort to available experimental values from human glioma models to ob- 
tain order-of-magnitude estimates of the intervening parameters in our equa- 
tions. First, the maximum cell density C*-*^^ has been estimated in previous 
works (see e.g. Rockne et al. 2010) to be of about 10^ cell/cm^. The oxygen 
concentration at which cells switch to hypoxic metabolism depends on the 
specific cell line but experimental evidences support for glioma the choice of 7 
mm Hg ( Vaupel|2004 1 . The Michaelis-Menten constant has to be smaller than 
this value y et larger than the anoxia threshold of about 0.7 mm Hg ( [Brown 



and Wilson 



et al.||2003 



20041. We have chosen it to be about O 



(T) _ 



Wilson|[2008| 



2.5 mm Hg (Dagu 



For the oxygen consumption parameters, although there is a high variabil- 
ity between different cell lines, we will take the values of U251 glioma cells 
from Griguer et al. (2008). As to the hypoxic oxygen consumption, this is a 
parameter difficult to assess. We will use data for the same cell line U251 but 
genetically modified to deplete its mitochondrial DNA so that the cell's oxida- 



tive capacity is virtually impaired. Data from Griguer et al. (20081 indicates 
that the normoxic cell uptake is five times higher than the hypoxic cell con- 
sumption. We have indeed assumed that the consumption rates were reported 
for normal conditions of pressure and temperature. Moreover, it is known that 



al. 



mm Hg (Kimura et 


Wilson 


2008 


) thus. 



60 



we will take the oxygen pressure at our boundary capillaries to be 
mm Hg. 

The oxygen diffusion coefficient is classically known to be around lO"*^ 
cm^/s ( Muelller-Klieser et al. 1984) while the cell diffusion coefficients are 



not so readily accessible to obtain in vivo. We resort to the data from |Waiig| 
et al. (2009), were the authors have exploited a recently developed tool for 



quantifying glioma net invasion rates in individual patients using routinely 
available magnetic resonance images (MRI) from 32 patients to derive a net 
measure of the invasive capacity of the overall dominant phenotype. Since 
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they obtain the parameters for the invasive phenotype, which is the most 
migratory one, and this is hnked to the hypoxic one in our model, we wih take 
their median migration rate 29 mm^/year as our fast diffusion coefficient for 
hypoxic cells. The normoxic cell diffusion will be taken to be smaller than the 
hypoxic one, in our case _D„ = O.lDh- 

Finally, the proliferation parameters can be estimated from typical dou- 
bling times that range from 24 to 48 h for glioma cells in vitro for well oxy- 
genated cells, while the mobile hypoxic cells doubling times are assumed to be 
longer (we have taken this parameter to range within two to four times the 
previous one). Table [l] summarizes the typical parameter values employed in 
our calculations. 



Table 1 Values of the biological parameters used into our 
Variable Description Value 



C(M) 



o. 



Maximum tumor 
cell density 



10" cell/cm^ 



2 Oxygen concentration 7 mmHg 

switch to hypoxia 



model equations. 
References 
i Rockne et al.|2010 l 

I Vaupel|2004 i 



2.5 mmHg ( Dagu et al.|2003 i 



O. 



Michaelis Menten 
oxygen threshold 



0.7 mmHg I Brown and Wilson|2004 i 



O. 



WT 



Oxygen level 
for anoxia 



02*1=1,2 Vessel oxigen 60 mmHg 
concentration 


Wilson|2008 


^.'l996j 


Kimura et a 




an Normoxic cell oxygen 10~^' Mol/(cs) I 


Griguer et al. 2008] 



consumption 



a„/5 Mol/(c s) I Griguer et al.|2008 l 



Hypoxic cell oxygen 
consumption 



| |Tiia and Moghe|2 002|~ 

c'vvang'c't^i: 'jxmv 



Diffusion coefficient for 
normoxic cells 



;.6-10-i"'cm^/s 



Diffusion coefficient for 10 Dn cm /s 





hypoxic cells 




( Wang et al. 2009 ) 


D02 


Oxygen diffusion 


10"'' cm^/s 


(Muelller-Klieser et al.|1984 





] Jewell and Gassmann|2001 1 



T„fi Phenotype switch time 1 h 

(normoxic to hypoxic) 

Thn Phenotype switch time 96 h 

(hypoxic to normoxic) 



Estimated 



Cell death mean 
time in anoxia 



48 h 



Estimated 



(Giesc et al."2003') 
(Berens and Giese 19991 
~ Giese et at 2003 1 * 



Normoxic cell doubling 
time 



24 h 



Hypoxic cell doubling 
time 



48 h 

0.1 Mol/mm^ 



Berens and Giese 1999|l 
Estimated 



AO2 



Sensitivity around 
transition threshold 



i Ma2umdar|1992[ | 



Oxygen exchange 
coefficient 



0.1 cm 
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2.4 Computational details 

The model given by Eqs. ([T]) and ^ describes the dynamics of cells in the bulk 
of the tumor under variable oxygenation conditions. As such, it can be studied 
on any particular geometry either regular or irregular with suitable boundary 
conditions and/or oxygen sources. It can also be coupled to detailed models of 



tumor vasculature such as those based on phase-field models (Travasso et al. 



2011) and/or combined with therapies having an oxygen-dependent effective- 



ness such as radiotherapy ( 


Jensen 


2009 


Vaupel 


2008 



Here, we wish to analyze our model to get insight on the formation of 
necrotic areas and hypercellular regions such as those observed in gliomas. To 
do so, we will stick to the simplest possible scenarios: either two dimensional 
rectangular tumor geometries with vessels located on two opposite boundaries 
and periodic boundary conditions on the other (opposed) two sides, or just one 
dimensional sections of the previous system with oxygen flow coming from the 
boundaries (end points) of a one-dimensional domain. 

The sources of nutrients and oxygen in the brain are the blood vessels 
which, in our 2D model, correspond to two of the four boundaries of the 
tumor tissue analyzed and, in the ID case, to two sides of the interval under 
study. Here, the blood vessels are assumed to be capillaries. Exchange of gases, 
nutrients, and waste materials occurs through the thin permeable walls of these 
capillaries. We will denote the oxygen concentration on the two lateral blood 
vessels at a given time t by (t) and (t). Oxygen flows from vessels to the 
tissue to balance the different oxygen pressures, thus the boundary conditions 
for the oxygen at the blood vessels are 

j^(0r(i)-02), * = 1,2. (4) 

The capillary length in our model simulations is taken to be < 1 mm to match 
the real capillary size ( Mazumdar|1992 ). In addition, oxygen concentration in 



the capillary network does not follow the variations induced by the circulation 
in the major blood vessels ( Dagu et al.|2003 ). Since our computational domain 



extends over a small spatial region (~700 fxm in length parallel to the blood 
vessels), we will assume the oxygenation to be the same independently on 
the position along the capillaries as our computational domain extends over a 
small spatial region (a strip of about 700 fim in length parallel to the blood 
vessels) . 

Likewise, the oxygen exchange coefficient Ji is considered to be constant 
and equal for both vessels, therefore | Ji| = | J2I = J- Bearing in mind that the 
average blood velocity for capillaries is around 0.1 cm/s ( [Mazumdar||1992 ) we 



can estimate its value to be J ^ 0.1 cm. Notice that our boundary conditions 
Q imply that there exists the possibility that oxygen spreads either away 
from or into the vessel should the tissue oxygen concentration be lower that 
or exceed the oxygen pressure at the vessel. 
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Once the computational domain has been fixed, we impose the bound- 
ary conditions for the cell densities to supplement model equations ([T]) and 
([2]) . The following three equations represent the flux conditions describing cell 
intravasation into the blood vessels 

(5a) 
(5b) 
(5c) 

Note that it is possible to calculate the number of cells that intravasate into 
the blood flow for t & [0,T], T being the studied temporal window, by means 
of the formula 

M{t) = I {JnCn + JhCh) dt, (6) 





1 


dx 


J n 


dCh 


1 


dx 


~ "a 


dCd 


= 0, 




dx 



where, in general, it is to be expected that the hypoxic intravasation coefficient 
Jh would be substantially larger than the normoxic one J„, due to the higher 
motility of the former cells. In metastatic processes, solid tumor cells enter 
into the lymphatic and blood vessels eventually leading to the formation of 



colonies in distant organs, thus Eqs. (5a) and (5b I represent the first stage of 



metastasis. However, GBM only rarely metastatize out of the brain. Although 
GBM cells may spread along blood vessel walls, primary brain tumors do 



neither intravasate to blood vessel walls (Bellail et al. 20041 nor invade into 
the bony calvarium which encloses the brain. Indeed, the capacity for vascular 
invasion is a primary requisite for wide hematogenous metastasis of malignant 
neoplasms. Thus, the absence of vascular invasion by primary brain tumors, 
taken together with the absence of lymphatics within the brain, may explain 
the observation that primary brain tumors only very rarely metastasize outside 
of the brain to ectopic organs. The unique patterns of invasion and metastasis 
shown by primary brain tumors suggest that these neoplasms are adapted to 
the microenvironment specifically present within the central nervous system. 

Therefore, we will assume henceforth that for gliomas Jn — Jh — 0. How- 
ever, our model equations together with Eqs. can be used in more general 
contexts, e.g. to elucidate the relation between acute hypoxia or cyclic hypoxia 
and metastatic processes in other types of tumors. 

To solve our model set of equations with the boundary conditions on the 
capillaries, plus periodic boundary conditions on the free sides, we have used 
standard finite differences (second order in space) either with an explicit fourth 
order Runge-Kutta method in time or an implicit second order in time Crank- 
Nicholson-type scheme. Typical parameter values for the spatial step were 
Ax — Ay— 5 /im, and At = 0.01 s for the explicit scheme. 

In all of the simulations to be presented later we thus assume an initial 
concentration of normoxic tumor cells in the region limited between two func- 
tional vessels. Thus, initially the hypoxic and necrotic cell densities are zero 
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for all X £ {0, Lx). The initial oxygen distribution is randomized with 24 mm 
Hg as the mean oxygen partial pressure as described by ( Vaupel et al.||1989 |. 
Except when otherwise stated, we will assume that the vessel Vi is disrupted at 
t = so that V2 remains as the only functional vessel through the simulation. 



3 Results 

3.1 Vaso-occlusive events trigger necrosis, the formation of pseudopalisades 
and vessel-cooption 

We have solved our model equations ([I]) and ([2| by including only the essential 
biological facts at the cellular level allowing us to understand the formation 
of palisading necrosis, and to get an insight on the sequence of events leading 
to their formation as well as estimates of their lifetimes and sizes. Once the 
parameters are chosen in the appropriate realistic ranges (see Table [T]) , the 
formation of pseudopalisades does not require a detailed fine-tuning of the 
parameters but it is a characteristic phenomenon arising in many of our simu- 
lations. It is interesting to highlight that pseudopalisades are transient states, 
being part of a sequence of events to be summarized in what follows. 

In Fig. [4] we present a first scenario of the results obtained when normoxic 
GBM cells are seeded close to one of two vessels located at the boundaries of 
a brain tissue slice of 300 ^m. Initial normoxic tumor cell density is maximum 
around one of the blood vessels (colored black in the lower side of Fig. |4]). The 
initial conditions for this simulations are taken to be of the form C„(x, 0) — 
^l-xy{Lj2)^ , X e (0,L,/2] and C„(a;,0) = for x e (ij2,L,] with 
Lj. —300 fim. We assume that the brain tissue is well oxygenated until the 
lower vessel is disrupted at t = so that, at the outset, the hypoxic and 
necrotic cell densities are zero in the whole computational domain. The results 
for the evolution of the different cell densities are summarized in Fig. |4] 

After the lower blood vessel undergoes a thrombus induced by the adja- 
cent high density of proliferating normoxic tumor cells, these cells experience 
an oxygen shortage and change their phenotype to hypoxic. This leads to an 
enhanced mobility and reduced proliferation in the perivascular region in an 
attempt to find a more nourishing environment to survive (Fig. [4]A). A com- 
plete phenotypic switch occurs in about 1 hour (see Fig. |4|3) so that the new 
hypoxic cells (Fig. |4]C) increase their flux to the upper vessel (in red) that is 
still functional, generating in the process a transient traveling wave of hypoxic 
cells (see Fig. |4|l) where the arrow shows the direction of the movement). 

Figs. |4]^E,F) depict the formation of a pesudopalisade in the scale of 1 day 
that develops during the subsequent days. A morphologic analysis of Figs. 
and F gives about 130 /xm as the characteristic width value. According to |Brat| 



et al. (20041, pseudopahsades with small widths have smaller necrotic cores. 
In our case, the necrotic core is about 30 /um (Fig. |4]3) and it is also formed 
in the same time scale of days. 
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Fig. 4 Spatio-temporal simulations showing the formation of necrotic areas sur- 
rounded by hypercellular regions (pseudopalisades). Initially, normoxic GBM cells 
are seeded close to the lower vessel. Pseudo color plots in A, B, C, E and G represent 
normoxic, normoxic (notice the scale), hypoxic, total (hypoxic +normoxic) populations and 
necrotic tissue densities, respectively. The horizontal and vertical axes correspond to time in 
days (except for B in hours) and space in fim, respectively. D and F depict the normoxic and 
total (normoxic + hypoxic) cell densities as a function of space for specific times indicated 
in each curve and also with vertical color lines in C and E. H shows the GBM cell number 
evolution with time for the different populations. At t = 0, the lower vessel suffers a throm- 
bus, while the upper vessel remains fully functional during the simulation. The parameters 
used in the simulations are listed in Table [T] 



Cell Waves around Necrotic Cores in GBM 



15 




vessel distance vessel distance 

Fig. 5 Initial conditions for the simulation displayed in Fig. [rj (A) Normoxic cell 
density. (B) Initial o xygen distribution w ith 24 mniHg as the average oxygen pressure for a 
healthy brain tissue ( [Vaupel et al.|1989| |. 




Fig. 6 Vascular damage as a driving force for migration in Glioblastoma Mul- 
tiforme. Formation of necrosis and the emergence of pseudopalisados is part of a series 
of related events that are crucial for the accelerated progression of high grade astrocytic 
gliomas. The development of an envelope of small scale pseudopalisados may play a relevant 
role in the global migratory dynamics and invasion of these malignant gliomas driven by 
vaso-occlusive events occurring at the small cellular scale, leading to the distinctive feature 
in GBM of showing significant necrotic areas where frequently central vessels are present 
l |Brat et al.|2004^ . 

Under persistent anoxic conditions, hypoxic cells located around the im- 
paired vessel start to die (Fig.|4]G;). Nevertheless, there is a significant density 
of migrating hypoxic cells that reach the functional vessel. Finally, Fig. |4fl 
illustrates the time evolution of the total GBM cell number for each popu- 
lation in the full simulated tissue. It is observed that hypoxic cell number 
increases significantly fast during the first hours after the thrombosis event 
but then moderates its growth rate. Normoxic cell numbers have a different 
growth pattern with an initial decay and a final relapse, once the functional 
vessel has been co-opted. So, in less than a week, tumor cells have invaded the 
300 fim of tissue. 

In a second set of simulations, displayed in Figure [7] we have studied the 
invasion of a tumor vessel by a pseudopalisading wave of hypoxic cells coming 
from another impaired vessel. As in the previous case, a high normoxic GBM 
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cell density is seeded around the lower blood vessel (see Fig. [7| , but now we 
assume that a small number of tumor (normoxic) cells have already migrated to 
the upper blood vessel starting to colonize the second vessel (400 fim distant 
from the lower vessel). The initial conditions for the tumor (normoxic) cell 
density and oxygen concentrations are shown in Fig. [5] and for the normoxic 
tumor cell density are given by C„(a:,0) = |(1 — tanh(^— y^)), x € {0,Lx/2] 
and Cnix,0) = 0.05(1 + tanh( "^~'^|~^" )) for x £ {L^/2,L^] with =400 ^m. 

The simulation starts when the lower vessel suffers a thrombus induced 
by the adjacent high density of proliferating normoxic tumor cells. The up- 
per vessel is fully functional during the simulation. This vessel is completely 
invaded, after 7 days, by a high tumor cell density that might induce new 
vaso-occlusive events (Fig. [Tj'V). The new hypoxic cells (Fig. ^p) generate a 
transient traveling wave (Fig. [7|l)) where the arrow shows the direction of the 
movement. Figures [7^ and F show the formation of a pseudopalisade in the 
time scale of one day that develops during the ensuing 7 days. 

A morphologic analysis of Figs. and F gives about 150 /im as the char- 
acteristic width value. The necrotic core is about 60 fim (Fig.jTjS) and it is also 
formed in the same time scale of days. Finally, it is observed in Fig. [7j3 that 
the hypoxic cell number increases significantly fast during the first hours after 
the thrombosis event but then moderates its growth rate. Normoxic cell num- 
bers have a different growth pattern with an initial decay and a final relapse, 
once the functional vessel has been co-opted. 

Qualitatively, it is important to underscore that this minimal model not 
only reproduces the clinical observations seen in immunohistochemical anal- 
ysis of histologic sections of GBM samples (Brat et al. 20041 but may also 
contribute to a better understanding of the origin and evolution of those struc- 
tures and to find unknown quantities such as the typical palisade lifetime or 
the relationship between palisade lifetime/size and vessel distance. Moreover, 
the macroscopic progression of GBM can be partially conceived, from a micro- 
scopic point of view, as consisting of a large number of intravascular throm- 
botic events (caused by unregulated tumor cell proliferation in the vicinity 
of the vessels) giving rise to the formation and coalescence of small necrotic 
foci coupled with subsequent episodes of hypoxic cell migrations in search of 
nearby functional blood vessels which will eventually suffer new intravascular 
thrombotic events. The envelope of many of these small-scale pseudopalisades 
contributes to the high density regions on larger spatial scales as those shown 
in Fig. [6) 



3.2 Palisading waves leads to invasion faster than that of pure random 
motion invasion 

In principle, it is not clear a priori how the vaso-occlusive events relate to 
the speed of the tumor invasion. Phenotypic characteristics of GBM such as 
hypercoagulation, hypoxia, and abnormal angiogenesis, may be linked at the 
molecular level, and hypoxia may coordinate the hypercoagulative activity of 
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Fig. 7 Spatio-temporal simulations showing the formation of perinecrotic pseu- 
dopalisades after a vaso occlusive event. In contrast with Figure |4] here both vessels 
are initially seeded with normoxic GBM cells. After 7 days, the upper vessel is completely 
invaded by a high GBM cell density and again a vaso occlusion event may occur at any 
moment. All of the subplots and the rest of parameters are as in Figure |4] and the initial 
data for oxygen and tumor cells are shown in Fig. [5] 
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Fig. 8 Palisading waves lead to invasion faster than that by pure random motion 
invasion. (A-B) Evolution of the tumor cell densities as a function of time and space 
when only the upper vessel is functional while the lower vessel becomes impaired at t = 0. 
Parameters are as in Fig. [4] Pseudocolor plot A depicts the total GBM cell density as a 
function of time, while plot B presents density curves at various times showing the advance 
of the tumor front. (C-D) Same as in plots (A-B) but when both vessels remain functional. 
All parameters and initial conditions are taken as in Figure |4] 



GBM cells and (protease-activated receptor) PAR-mediated angiogenic signal- 
ing dSvensson et al.||201lD . 

On the one hand, one might naively argue that the vessel maintenance 
would lead to an increased number of cells that by random motion (diffusion) 
would give rise to a faster invasion due to a higher tumor cell number. On the 
other hand, one might also argue that vaso-occlusive events, despite result- 
ing in hypoxia, cell death and slower proliferation, could induce more mobile 
cells and might instead accelerate the overall GBM progression. To gain fur- 
ther quantitative insight of what is the real situation, i.e. what is the role of 
vaso-occlusive events in relation with the tumor progression speeds, we have 
performed a series of numerical simulations. 

Fig.|8]is a representative example of our results. Subplots A,B in Fig.[8]cor- 
respond to the case discussed previously where the upper vessel works normally 
but the lower vessel becomes impaired due to the formation of a thrombus. 
On the contrary, Fig. |8p,D shows the tumor progression when both vessels 
remain functional. In the first case, the functional vessel is invaded by the tu- 
mor after 8 days by a high tumor cell density (more than 45%) however when 
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no thrombus is formed and both vessels stay functional (subplots C and D) 
invasion becomes significantly slower; requiring 13.5 days to attain the same 
tumor cell density around the upper vessel. Fig. |8] displays the coexistence of 
two mechanisms driving invasion: the first one being purely diffusive and simi- 
lar to that governed by a free growth (when resources are not limited) obeying 
the Fisher-Kolmogorov equation, as described by (Swanson et al. 2008) and 
related works. The second one shows an accelerated invasion caused by micro 
thrombi formations in the tumor vessels. Thus, on the basis of our model, it 
appears that vaso-occlusive events in GBM may play a key role in accelerating 
the tumor invasion process through pseudopalisade formation. 

Thus our model simulations suggest that the use of chemical agents slowing- 
down the vessel impair might delay tumor progression in GBM patients. This 
fact will be discussed throughly in Sec. |4j 



3.3 There is a relation between the size of broken vessels, the palisade 
lifetime and the width of the necrotic region 

To further substantiate our results, we have carried out a numerical study of 
the relation between the separation between a thrombotic vessel and another 
functional vessel, the palisade width and its lifetime, as depicted in Fig. [9] The 
characteristic width of the pseudopalisade corresponds to the dimensions of 
the hypercellular region and its lifetime is the time elapsed since the palisade 
is formed until it disappears. Although the standard distances between blood 
vessels are in the range of 100 to 400 ^m, due to the irregular and aberrant 
structure of the tumor vasculature, we have considered a wider range (from 50 
to 700 fj,m) assuming that those histologic samples showing distances larger 
than 700 fim will be formed by cells for which their main oxygen source would 
probably come from a vessel close to them but not visible in pathological 
analyses due to its spatial arrangement. 

Although our analysis correspond to one-dimensional scenarios we still 
expect to get a hint on the main tendencies, since more complex higher- 
dimensional scenarios with a realistic three-dimensional vessel network are 
too complicated to be sistematized. 

In our computational window, the distance between blood vessels equals 
the tissue size and so when this is smaller than 100 fim the formation of 
pseudopalisades is practically non-existent. In addition. Fig. [9] shows how the 
pseudopalisade size depends on the distance between blood vessels (that may 
be also expected to be roughly proportional to the vessel size), growing in a 
linear fashion as the distance increases until it reaches a plateau value close 
to 200 /im. Notice that vessel separations larger than 500 /im give rise to 
pseudopalisades that do not experience a significant width change despite the 
fact that a smaller difference between the pseudopalisade cell density and the 
adjacent tissue cell density is observed in our calculations. Also, as it might be 
expected, there is no appreciable necrosis until the distance between vessels is 
larger than 150 /im. Then, in the range from 150 to 450 /im the necrotic core 
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Palisade size, lifetime and necrotic core 




Blood vessels distance (^m) 

Fig. 9 Dependences between thrombotic-functional vessel separation, palisade 
width and its lifetime. In all of our simulations GBM normoxic cells are seeded close 
to the impaired vessel with a small fraction of them already invading the functional vessel. 
Initial conditions are chosen similar to those of Fig.[7]and proportional to the tissue size. At 
t = one of the vessels becomes thrombotic while the second vessel remains fully functional 
during the simulated time window. Characteristic pseudopalisade width and lifetime data 
are labelled by blue triangles, and red squares, respectively. Also shown is the dependence 
of the necrotic core width with the tissue size (green circles). The horizontal axis represents 
the distance between blood vessels (tissue size) in fim. The vertical axes label (left one) the 
palisade and the necrotic foci widths in /im; (right one) the pseudopalisade lifetime in hours. 
The parameters used in the simulations are listed in Table [T] 



grows exponentially and from here its growth becomes linear (Fig. |9| . In other 
words, if the distance from a thrombotic vessel to the nearest functional vessel 
is short enough, the cells arrive to better oxygenated areas without experienc- 
ing severe hypoxic stresses (thus diminishing their death) and reducing the 
palisade lifetime. However, if this distance is larger than 150-200 /im (or if the 
broken vessels are larger), there will be a considerable number of hypoxic cells 
unable to scape from severe hypoxia giving rise to larger necrotic areas and 
longer palisade lifetimes. Here again, the model although biologically simple 
and solved on a simple domain, allows to explain the observed feature that 
large pseudopalisades, i.e. pseudopalisades surrounding larger broken vessels 
are more persistent and thus easier to observe, than small ones (cf. Fig. ph. 
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4 Perspectives and Therapeutic implications 

4.1 Perspectives 

The model presented in this work can be extended in many ways. However 
the addition of vasculature seems the more natural way providing the basis 
to construct a more complex in silico framework to apply the different thera- 
peutic modalities in a more realistic manner once the underlying mechanisms 
involved in palisade formation have been laid down. One possibility is to main- 
tain equations (la), (lb I and (Ic) respectively for the normoxic, hypoxic and 



necrotic populations substituting the microscopic oxygen densities by macro- 
scopic variables such as blood vessel densities or vessel functionality v{x, t) 
that are to be related to the oxygenation conditions. This vessel density would 
evolve according to an equation of the form 



dv 
dt 



1 



1 

Tt 



THROMBOSIS (C„). 



(7) 



accounting for the fact that, given that space is available, hypoxic-migratory 
cells C/i contribute to the generation of the new vasculature v through the 
secretion of vascular endothelial grow factors (VEGF) (first term in Eq. ^ 
and that normoxic cells contribute to the vaso-occlusive events (second term 
in Eq.[7]). 

These ideas might provide a simple way to incorporate palisade forma- 
tion due to micro-vascular damage to full-brain scale mathematical models 



of glioma progression ( 


Frieboes et al. 


2007 


Swanson et al. 


2008 


Bondiau et 


al. |2008| 


Eikenberry et al. 2009| Konukoglu et al. 2010 Rockne et al. 2010 


Perez-Garci'a et al.|2011) may allow to construct more realistic "macroscopic" 



models. These models might be better suited to describe more realistically 
glioma progression and incorporate the effect of therapies and as such to pro- 
vide clinicians with valuable tools for treatment optimization. 

Other possible extensions of our model regard the presence of a small 
number of treatment-resistance (cancer-stem) cells in glioblastoma multiforme 
(GBM) and whether their metabolism can be effectively targeted. Indeed, de- 
termining whether this small fraction mainly depends on aerobic glycolysis 
or on oxidative phosphorylation has recently been the subject of considerable 
interest (PistoUato et al. 2010). It is known that in glycolytic tumors, which 



is the case of high grade gliomas, oxidative phosphorylation is not completely 
shut down since it provides a significant amount of ATP even under low glu- 
cose requirements ( Vander Heiden et al.||2009 Mathupala et al.|[2010 |. In this 
work we have focused on oxygen as one key factor known to strongly affect 
GBM progression and response to treatment (e.g. radiotherapy). While our 
distinction of the tumor cell phenotypes has relied on their oxic state, we an- 
ticipate that, among the migratory hypoxic cells seen in the pseudopalisades, 
a small fraction of them is formed by treatment-resistance cells. Actually, our 
simulations suggest that during the time window where hypoxic cells migrate 
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from an impaired blood vessel to a functional one any treatment critically de- 
pending on the presence of oxygen will exhibit a poor clinical response. The 
importance of oxygen in these subpopulations is further supported by recent 
experiments by Vlashi et al. (20111 in which it was observed that in various 
GBM cell lines both progenitor cells as well as those exhibiting a stem-cell- 
like behavior consume less glucose and produce less lactate while maintaining 
higher ATP levels than their differentiated progeny. Their conclusions were 
that GBM stem cells rely mostly on oxidative phosphorylation. Even from an 
evolutionary point of view, the collective migration of hypoxic GBM cells when 
forming pseudopalisades appears to be a selection mechanism of the fittest cell 
phenotypes. 

It is necessary to underscore that perinecrotic pseudopalisades formation 
in GBM is not a minor detail or yet another layer to add to the complexity 
of high grade glioma growth. There is ample evidence that the formation of 
necrosis and the appearance of pseudopalisades is part of a series of related 
events that are crucial for the accelerated progression of high grade astrocytic 
gliomas. Thus, the development of an envelope of small scale pseudopalisades 
may play a relevant role in the global migratory dynamics and invasion of 
these malignant gliomas driven by vaso-occlusive events occurring at the small 
cellular scale, leading to the distinctive GBM feature of showing significant 
necrotic areas, as observed in clinical imaging as it is exhibited in Figure [6j 

On the one hand, this hystopatologic feature is a distinctive feature of 
high-grade gliomas and does not arise in other types of tumors, probably due 
to the very high mobility of tumor cells in astrocytoma and the low resistance 
to invasion by the brain parenchyma. On the other hand, it is well known that 
gliomas are very complex and heterogeneous and yet, they can be understood 
in terms of two cell subpopulations corresponding to two dominant pheno- 
types ([DeBerardinis et al.||2007[ |Giese et al!]|2003| [Keunen et"al]|2011[ |Onishi 



et al.|[20TT ): a population of proliferative, less mobile cells located in perivas 



cular areas and a second population of less proliferative and more mobile cells 
typically infiltrating the brain parenchyma. This means that minimal math- 
ematical models such as those based on Fisher-Kolmogorov type equations 



accounting for GBM progression (Frieboes et al. 2007 


Swanson et al. 2008 


Bondiau et al. 


2008 


Eikenberry et al. 2009| Konukog 


lu et al. 2010 


Rockne 


et al.pOlO 


Perez-Garci'a et al.||2011 


I may benefit from incorporating the two 



different phenotypes and necrosis as it is done in this paper. 

However, in order to incorporate necrosis not only phenomenologically, as 
in Perez-Garci'a et al. (2011), but including at least some underlying mecha- 



nisms, it is necessary to take into account the basic intervening steps involved 
in palisade formation. In fact, to allow for the vasculature growth in response 
to the pro-angiogenic factors and breakup due to vaso-occlusive events, blood 
vessel dynamics must be included in the model. Moreover, oxygen variations 
occur in very small spatial scales and at very fast time scales that cannot be 
accounted for in detail when solving the problem at larger full-brain scales and 
for times of the order of months. This means that more global models must 
rely on the higher scale effective vascular density or vascular capacity m{x,t) 
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rather than on oxygen. This is essential since one of the conclusions of Brat et 



al. (2004) and Brat and Van Meir (2004), also supported by our simulations, 



is that vascular damage is not only a key event in the formation of necrotic 
areas but also a driving force for migration, and this dynamics is linked to the 
evolution of the cell densities that are responsible for the vessel destruction 
(normoxic cell population in the perivascular areas) and creation (hypoxic cell 
subpopulations secreting VEGF and other pro-angiogenic factors). 



4.2 Therapeutic implications 



One of the main results of our work is the prediction that targetting vaso- 
occlusive events in GBM might delay the tumor progression, a result with 
relevant potential implications since GBM is a type of tumor for which even 
modest survival increases are considered to be a big success because of its poor 
prognosis. 

Thromboembolism, i.e. vaso-occlusive events of large vessels, is well rec- 
ognized as a major complication of cancer and known to be a common cause 
of death in cancer patients. There is strong evidence linking venous throm- 
boembolic events and malignancy. Laboratory markers of coagulation activa- 
tion such as thrombinantithrombin complex or prothrombin fragments 1 + 
2 support the premise that malignancy is a hypercoagulable state. Inflamma- 
tory cytokines (e.g. tumor necrosis factor and interferon), coagulation proteins 
(e.g. tissue factor and factor VIII), and procoagulant microparticles may be 
elevated in patients with malignancy. However, the relative contribution of 
chemotherapeutics, tumor cells, endothelium, and circulating procoagulants 



in promoting thrombus formation continues to be investigated (Khorana and 



Francis] [2008{ |Furie and Fmie| [20061 [Zwicker et ar][2007| [Green and Kwaan 
2009| 



Thromboembolic complications include a broad spectrum of clinical prob- 
lems, a fact that has lead to the use of thromboprophylaxis in a variety of 



forms (Dolovich et al. 2000[ Khorana and Francis 2008 Green and Kwaan 



2009 ) for cancer patients. Specifically, gliomas have a high incidence of venous 
thromboembohsm (VTE); with several studies suggesting that 25%-30% of 
these patients sustain thromboembolic events ( Streiff et al.[[2004 Simanek et 



al. 20071. The fact that glioma cell lines secrete pro-thrombotic factors has 



been known for a long time (Bastida et al. 1984) and it is also known that 



the more tumoral tissue is removed during surgery of high-grade gliomas the 



less- likely are the patients to die from VTE (Brose and Lee 2008 Simanek 
et al.[[2007 ). This fact has led to the consideration of thromboprophylaxis for 
glioma patients, specially in subjects with higher potential risks of developing 



venous thromboses (Hamilton et al.[[1994 Batchelor and Byrne[[2006 


et al.[[2010 


Khorana and Francis[[2008 


). 



In addition to the prevention of VTE there are several direct mechanisms 
of action of some anticoagulants, such as LMWH, on tumor cells: direct cell 
killing ( Santos et al.[[2010 l, anti- angiogenic effects ( Svensson et al.[[2011 ) and 
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many other (see e.g. Green and Kwaan 2009 Chap. 15). The results from 



our model seem to imply that an additional indirect antitumoral effect from 
thromboprophyylaxis might be expected related to the delay of tumor invasion, 
different from the direct antitumoral effect or the increase of life expectancy 
obtained from the prevention of the formation of big coagulates. 

The main limitation in using LMWH in post-operative patients has been 
related to the possible increase in bleeding events, that in the case of brain 
surgery are potentially very harmful. However, a limited phase II clinical study 
of LMWH therapies for high grade glioma patients by Robins et al. ( 2008 ) has 
shown that the therapy is safe. A recent phase HI trial with a small number 
of patients seems to confirm this trend Perry et al. (2010). 

If, as predicted, preventing capillar thrombosis induced by the tumor, re- 
sults in a growth delay for the tumor, a synergistic positive effect with radiation 
therapy is to be expected. Having a better oxygenation and reducing hypoxia 
would lead to an enhanced tumor radio-sensitivity, that might provide an extra 
benefit in survival. 

Whether these predictions are clinically relevant remains to be studied, first 
in animal models testing if the use of LMWH induces a significant reduction of 
the formation time for palisades and delays tumor progression and, in case of 
success, in more complete clinical trials. The fact that mathematical modelling 
might contribute to designing therapies against such an aggressive type of brain 
tumor is encouraging and we expect our work to stimulate more theoretical 
and experimental studies to follow in this direction. 



5 Conclusions 

In this work we have put forward a simple continuous biomathematical model 
allowing for the description of the formation of hypoxic hypercellular regions 
around necrotic cores in glioblastoma multiforme: the so-called pseudopal- 
isades. Pseudopalisading necrosis and microvascular hyperplasia are two of 
the most powerful predictors of poor prognosis among invasive gliomas. Rather 
than mere morphologic markers, these structures are almost certainly mech- 
anistically instrumental in the accelerated growth properties that character- 
ize the transition from a lower to a higher grade astrocytic glioma, such as 
glioblastoma. 

Our biomathematical model quantifies the migration process influenced by 
the phenotypic switch of glioblastoma cells under hypoxic conditions caused 
by vaso-occlusive episodes. In spite of the fact that tumor forms a three di- 
mensional structure with a very complex dynamics, our numerical simulations 
clearly evidence that it is possible to observe the dependence between palisade 
lifetime, its necrotic core and characteristic width related to the distance from 
a thrombotic to a functional vessel within a minimal mathematical model. 

We have also provided quantitative metrics of the characteristic times of 
formation and lifetimes of pseudopalisades after a vaso-occlusion event as well 
as their width and necrosis. We stress that these important metrics are very 
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challenging to measure in vivo because of the small spatial scale of the phenom- 
ena involved and the fact that biopsies do not provide dynamical information 
in time but only a fixed image for a given time. In addition, our results in- 
dicate that smaller pseudopalisades would show shorter lifetimes and reduced 
necrotic cores while larger pseudopalisades would display longer lifetimes until 
reaching a threshold around the 8th day after the vaso-occlusion event. From 
our calculations, pseudopalisades tend to present a characteristic width of 200 
fiin when blood vessel distances are larger than 500 fiui (within the range 
explored up to 700 fim). 

Furthermore, our simulations show that palisading waves lead to a faster 
invasion than that by a pure random motion invasion with unlimited resources, 
implying that the mechanism of vessel break up might accelerate the glioblas- 
toma progression. This suggests that targeting vaso occlusion events, probably 
by means of anticoagulants such as low molecular weight heparine, or other 
antithrombotic agents, might be of potential interest for therapeutics aimed 
at delaying glioblastoma growth. 
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